# table a7 left, appendix 

wru_freqs<-read.csv( "wru_freq.csv")
wru_freqs<- aggregate(matched~wrurace, data=wru_freqs, FUN=mean)
xtable(wru_freqs, digits=3)


#how many white people do we have to switch from 1 -> 0 to get no stat sig?  (table a7 right)
set.seed(7654321)
results_df<- data.frame(0:50)
results_df$reduc_white<-NA
results_df$t_val<-NA

css2<- read.csv("css05_cv.csv")
css2$part<- css2$cash + css2$vouchers
css2$part[css2$part==2]<-1
css2$part<- css2$part*100
css2$voucher_yr<-0
css2$voucher_yr[css2$elec_year>2016]<-1
css2$iswhiteperson<- 0
css2$iswhiteperson[css2$race=="European"]<- 1
css2$iswhiteperson[is.na(css2$race)]<- NA
css2$iswhiteperson[css2$race==""]<- NA

css2$income<- as.numeric(substr(css2$income,2,length(css2$income)))
css2$income_perc<- ceiling(ecdf(css2$income)(css2$income)*10)
css2$highinc<-0
css2$highinc[css2$income>100000]<-1
css2$highinc[is.na(css2$income)]<-NA
css2$highvote<-0
css2$highvote[css2$sums>13]<-1
css2$isdem<- 0
css2$isdem[css2$party=="Democratic"]<- 1
css2$isdem[is.na(css2$party)]<- NA



qqq1<- summary(felm(part~iswhiteperson*voucher_yr|as.factor(id)+as.factor(elec_year)|0|0 , data=css2),robust=T)
i=1
results_df$reduc_white[i]<- 0
results_df$t_val[i]<- qqq1$coefficients[3,3]
results_df$coef<-NA
results_df$coef[i]<-qqq1$coefficients[3,1]

css3a<- css2[css2$iswhiteperson!=1 | css2$part==0,]
css3b<- css2[css2$iswhiteperson==1 & css2$part==100,]
for (i in 2:25) {
  r<- i-1 
  binom_prob <- 1- (r *.02)
  results_df$reduc_white[i]<- binom_prob
  css3b$part<- rbinom(47616, 1,binom_prob) *100
  css3<- bind_rows(css3a, css3b)
  qqq1<- summary(felm(part~iswhiteperson*voucher_yr|as.factor(id)+as.factor(elec_year)|0|0 , data=css3),robust=T)
  results_df$t_val[i]<- qqq1$coefficients[3,3]
  results_df$coef[i]<-qqq1$coefficients[3,1]
  print(i)
}
print(xtable(results_df[1:15,2:4], digits=3), include.rownames=FALSE)




#######################################################################################
#appendix table 6
#####################################################################################

set.seed(7654321)
results_df_white<- data.frame(0:40)
results_df_white$t_val<-NA
results_df_white$coef_white<-NA


qqq1<- summary(felm(part~iswhiteperson*voucher_yr|as.factor(id)+as.factor(elec_year)|0|0 , data=css2),robust=T)
i=1
results_df_white$reduc_v[i]<- 1
results_df_white$t_val[i]<- qqq1$coefficients[3,3]
results_df_white$coef_white[i]<- qqq1$coefficients[3,1]

css3a<- css2[css2$vouchers==0,]
css3b<- css2[css2$vouchers==1,]
for (i in 2:20) {
  binom_prob <- 1- (i *.01)+.01
  results_df_white$reduc_v[i]<- binom_prob
  css3b$vouchers<- rbinom(NROW(css3b), 1,binom_prob)
  css3<- bind_rows(css3a, css3b)
  css3$part<- css3$cash + css3$vouchers
  css3$part[css3$part==2]<-1
  css3$part<- css3$part*100
  qqq1<- summary(felm(part~iswhiteperson*voucher_yr|as.factor(id)+as.factor(elec_year)|0|0 , data=css3),robust=T)
  results_df_white$t_val[i]<- qqq1$coefficients[3,3]
  results_df_white$coef_white[i]<- qqq1$coefficients[3,1]
  print(i)
}


results_df_white$income<-NA
results_df_white$coef_income<-NA
css3a<- css2[css2$vouchers==0,]
css3b<- css2[css2$vouchers==1,]
set.seed(1234567)
for (i in 1:20) {
  binom_prob <- 1- (i *.01)+.01
  css3b$vouchers<- rbinom(NROW(css3b), 1,binom_prob)
  css3<- bind_rows(css3a, css3b)
  css3$part<- css3$cash + css3$vouchers
  css3$part[css3$part==2]<-1
  css3$part<- css3$part*100
  qqq1<- summary(felm(part~income_perc*voucher_yr|as.factor(id)+as.factor(elec_year)|0|0 , data=css3),robust=T)
  results_df_white$income[i]<- qqq1$coefficients[3,3]
  results_df_white$coef_income[i]<- qqq1$coefficients[3,1]
  print(i)
}

results_df_white$voting<-NA
results_df_white$coef_voting<-NA
set.seed(2468)
for (i in 1:20) {
  binom_prob <- 1- (i *.01)+.01
  css3b$vouchers<- rbinom(NROW(css3b), 1,binom_prob)
  css3<- bind_rows(css3a, css3b)
  css3$part<- css3$cash + css3$vouchers
  css3$part[css3$part==2]<-1
  css3$part<- css3$part*100
  qqq1<- summary(felm(part~sums*voucher_yr|as.factor(id)+as.factor(elec_year)|0|0 , data=css3),robust=T)
  results_df_white$voting[i]<- qqq1$coefficients[3,3]
  results_df_white$coef_voting[i]<- qqq1$coefficients[3,1]
  print(i)
}

results_df_white$gop<-NA
results_df_white$coef_gop<-NA
css2$republican<- 0
css2$republican[css2$party=="Democratic"]<- 0
css2$republican[css2$party=="Republican"]<- 1
css2$republican[css2$party=="Non-Partisan"]<- 0
css3a<- css2[css2$vouchers==0,]
css3b<- css2[css2$vouchers==1,]
set.seed(8642)
for (i in 1:20) {
  binom_prob <- 1- (i *.01)+.01
  css3b$vouchers<- rbinom(NROW(css3b), 1,binom_prob)
  css3<- bind_rows(css3a, css3b)
  css3$part<- css3$cash + css3$vouchers
  css3$part[css3$part==2]<-1
  css3$part<- css3$part*100
  qqq1<- summary(felm(part~republican*voucher_yr|as.factor(id)+as.factor(elec_year)|0|0 , data=css3),robust=T)
  results_df_white$gop[i]<- qqq1$coefficients[3,3]
  results_df_white$coef_gop[i]<- qqq1$coefficients[3,1]
  print(i)
}

out_tab<- results_df_white[,c("reduc_v","coef_white","t_val","coef_income","income","coef_voting","voting","coef_gop","gop")]
print(xtable(out_tab[1:20,], digits=3), include.rownames=FALSE)



#####################################################################
############################  fig 1 (main text)
##
##
##
#####################################################################

css2<-read.csv("css11_cv.csv")

css2<- subset(css2, css2$elec_year>2015)


polpart_tab<-aggregate(vouchers~sums, data=css2, mean)

age_tab<-aggregate(vouchers~age_at_elecyr, data=css2, mean)

css2$estHHinc<- as.numeric(paste0(substr(css2$income,2,nchar(css2$income))))
css2$income_perc<-ceiling(ecdf(css2$estHHinc)(css2$estHHinc)*10)
income_tab<-aggregate(vouchers~income_perc, data=css2, mean)



ethn_tab<-aggregate(vouchers~race, data=css2, mean)

party_tab<-aggregate(vouchers~party, data=css2, mean)



gg_party<- ggplot(party_tab, aes(x= factor(party), y=vouchers*100))+geom_point(size=2)+geom_line()+
  xlab("")+ylab("Percent using\n Democracy Vouchers")+
  ggtitle("Party")+theme_light()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none", text = element_text(size = 20)) 
#gg_party



gg_inc<- ggplot(income_tab, aes(x= income_perc, y=vouchers*100))+geom_point(size=2)+geom_line()+
  xlab("")+ylab("Percent using\n Democracy Vouchers")+
  ggtitle("Income Decile")+
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+theme_light()+ theme(text = element_text(size = 20)) 
#gg_inc


age_tab$age<- as.numeric(paste0(age_tab$age_at_elecyr))
age_tab<- subset(age_tab, age_tab$age<95)
age_tab<- subset(age_tab, age_tab$age>29)
gg_age<- ggplot(age_tab, aes(x= age, y=vouchers*100))+geom_point(size=2)+geom_line()+
  xlab("")+ylab("")+
  ggtitle("Age")+theme_light()+ theme(text = element_text(size = 20)) 
gg_age

ethn_tab<-ethn_tab[2:5,]
ethn_tab$ethn<- c("Asian","White","Hispanic","Black")




hpp<- subset(css2, css2$sums>12)
inc_polpart<-aggregate(vouchers~income_perc, data=hpp, mean)
gg_incpolpart<- ggplot(inc_polpart, aes(x= income_perc, y=vouchers*100))+geom_point(size=2)+geom_line()+
  xlab("")+ylab("Percent using\n Democracy Vouchers")+
  ggtitle("Income Decile for Frequent Voters")+
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10))+theme_light()+ theme(text = element_text(size = 20))   
gg_incpolpart

ethn_tab2<-aggregate(vouchers~race, data=css2[css2$sums>12,], mean)
ethn_tab2<-ethn_tab2[2:5,]
ethn_tab2$ethn<- c("Asian","White","Hispanic","Black")
gg_eth2<- ggplot(ethn_tab2, aes(x= factor(ethn), y=vouchers*100))+geom_point(size=2)+geom_line()+
  xlab("")+ylab("")+
  ggtitle("Ethnorace for Frequent Voters")+theme_light()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none",text = element_text(size = 20)) 


party_tab2<-aggregate(vouchers~party, data=css2[css2$sums>12,], mean)
gg_party2<- ggplot(party_tab2, aes(x= factor(party), y=vouchers*100))+geom_point(size=2)+geom_line()+
  xlab("")+ylab("")+
  ggtitle("Party for Frequent Voters")+theme_light()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none", text = element_text(size = 20)) 





legend_b <- get_legend(
  gg_age +
    guides(color = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom")
)


gg_eth<- ggplot(ethn_tab, aes(x= factor(ethn), y=vouchers*100))+geom_point(size=2)+geom_line()+
  xlab("")+ylab("Percent using\n Democracy Vouchers")+
  ggtitle("Ethnorace")+theme_light()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none",text = element_text(size = 20)) 


gg_polpart<- ggplot(polpart_tab, aes(x= sums, y=vouchers*100))+geom_point(size=2)+geom_line()+
  xlab("")+ylab("Percent using\n Democracy Vouchers")+
  ggtitle("Past Elections Voted")+
  scale_x_continuous(breaks = scales::pretty_breaks(n = 14))+theme_light()+ theme(text = element_text(size = 20)) 


prow2<-plot_grid(gg_polpart,gg_age, gg_inc, gg_incpolpart,  gg_party, gg_party2, gg_eth,gg_eth2, ncol=2)
prow2



